
clear

set matsize 2000

* ssc install runmlwin

use TimelinePartiesMI_New.dta, clear

* RENAME VARIABLES TO MATCH SCRIPT BELOW
drop daysbeforeED
rename daysbeforeED_0_001 daysbeforeED

drop t_*
rename tk_* t_*

keep log_ae daysbeforeED pr t_pr partycentric t_partycentric enpp t_enpp vote_ t_vote partyold t_partyold niche t_niche gov_ t_gov inc t_inc cipollse countryid partyid _*_log_ae

gen cons=1
gen re1=1
gen re2=1
gen id=_n

order log_ae daysbeforeED pr t_pr enpp t_enpp partycentric t_partycentric vote_ t_vote partyold t_partyold niche t_niche gov_ t_gov inc_ t_inc cipollse countryid partyid cons re* id _*_log_ae

sort countryid partyid

save TimelinePartiesMLWin.dta, replace

*** MULTIPLE IMPUTATION OF MLWin ESTIMATES ***

global MLwiN_path "C:\Program Files\MLwiN v3.04\mlwin.exe"

*matrix list e(b)
*matrix list e(V)

gen M=50
gen m=.

gen b_daysbeforeED=.
gen b_enpp=.
gen b_t_enpp=.
gen b_t_partycentric=.
gen b_vote=.
gen b_t_vote=. 
gen b_partyold=. 
gen b_t_partyold=. 
gen b_t_niche=. 
gen b_gov=. 
gen b_t_gov=. 
gen b_inc=.
gen b_t_inc=. 
gen b_cipollse=. 
gen b_cons=.
gen b_level1=.
gen b_level2=.
gen b_level3=.
gen b_parl=.
gen b_t_parl=.
gen b_pr=.
gen b_t_pr=.

gen se_daysbeforeED=.
gen se_enpp=.
gen se_t_enpp=.
gen se_t_partycentric=.
gen se_vote=.
gen se_t_vote=. 
gen se_partyold=. 
gen se_t_partyold=. 
gen se_t_niche=. 
gen se_gov=. 
gen se_t_gov=. 
gen se_inc=.
gen se_t_inc=. 
gen se_cipollse=. 
gen se_cons=.
gen se_level1=.
gen se_level2=.
gen se_level3=.
gen se_parl_=.
gen se_t_parl=. 
gen se_pr_=.
gen se_t_pr=. 

gen V_daysbeforeED=.
gen V_enpp=.
gen V_t_enpp=.
gen V_t_partycentric=.
gen V_vote=.
gen V_t_vote=. 
gen V_partyold=. 
gen V_t_partyold=. 
gen V_t_niche=. 
gen V_gov=. 
gen V_t_gov=. 
gen V_inc=.
gen V_t_inc=. 
gen V_cipollse=. 
gen V_cons=.
gen V_level1=.
gen V_level2=.
gen V_level3=.
gen V_parl=.
gen V_t_parl=.
gen V_pr=.
gen V_t_pr=.

rename vote_ vote
rename gov_ gov
rename inc_ inc

foreach x of numlist 1/50 {

replace m=`x' in `x'

runmlwin _`x'_log_ae daysbeforeED cipollse pr t_pr enpp t_enpp vote t_vote partyold t_partyold t_niche gov t_gov inc t_inc cons, level3(countryid: cons) level2(partyid: cons) level1(id: cons) nopause

matrix bts = e(b)
scalar bt = bts[1,1]
replace b_daysbeforeED=bt in `x'
scalar bt = bts[1,2]
replace b_cipollse=bt in `x'
scalar bt = bts[1,3]
replace b_pr=bt in `x'
scalar bt = bts[1,4]
replace b_t_pr=bt in `x'
scalar bt = bts[1,5]
replace b_enpp=bt in `x'
scalar bt = bts[1,6]
replace b_t_enpp=bt in `x'
scalar bt = bts[1,7]
replace b_vote=bt in `x'
scalar bt = bts[1,8]
replace b_t_vote=bt in `x'
scalar bt = bts[1,9]
replace b_partyold=bt in `x'
scalar bt = bts[1,10]
replace b_t_partyold=bt in `x'
scalar bt = bts[1,11]
replace b_t_niche=bt in `x'
scalar bt = bts[1,12]
replace b_gov=bt in `x'
scalar bt = bts[1,13]
replace b_t_gov=bt in `x'
scalar bt = bts[1,14]
replace b_inc=bt in `x'
scalar bt = bts[1,15]
replace b_t_inc=bt in `x'
scalar bt = bts[1,16]
replace b_cons=bt in `x'
scalar bt = bts[1,17]
replace b_level3=bt in `x'
scalar bt = bts[1,18]
replace b_level2=bt in `x'
scalar bt = bts[1,19]
replace b_level1=bt in `x'

matrix bts = e(V)
scalar bt = bts[1,1]
replace se_daysbeforeED=sqrt(bt) in `x'
scalar bt = bts[2,2]
replace se_cipollse=sqrt(bt) in `x'
scalar bt = bts[3,3]
replace se_pr=sqrt(bt) in `x'
scalar bt = bts[4,4]
replace se_t_pr=sqrt(bt) in `x'
scalar bt = bts[5,5]
replace se_enpp=sqrt(bt) in `x'
scalar bt = bts[6,6]
replace se_t_enpp=sqrt(bt) in `x'
scalar bt = bts[7,7]
replace se_vote=sqrt(bt) in `x'
scalar bt = bts[8,8]
replace se_t_vote=sqrt(bt) in `x'
scalar bt = bts[9,9]
replace se_partyold=sqrt(bt) in `x'
scalar bt = bts[10,10]
replace se_t_partyold=sqrt(bt) in `x'
scalar bt = bts[11,11]
replace se_t_niche=sqrt(bt) in `x'
scalar bt = bts[12,12]
replace se_gov=sqrt(bt) in `x'
scalar bt = bts[13,13]
replace se_t_gov=sqrt(bt) in `x'
scalar bt = bts[14,14]
replace se_inc=sqrt(bt) in `x'
scalar bt = bts[15,15]
replace se_t_inc=sqrt(bt) in `x'
scalar bt = bts[16,16]
replace se_cons=sqrt(bt) in `x'
scalar bt = bts[17,17]
replace se_level3=sqrt(bt) in `x'
scalar bt = bts[18,18]
replace se_level2=sqrt(bt) in `x'
scalar bt = bts[19,19]
replace se_level1=sqrt(bt) in `x'

matrix bts = e(V)
scalar bt = bts[1,1]
replace V_daysbeforeED=bt in `x'
scalar bt = bts[2,2]
replace V_cipollse=bt in `x'
scalar bt = bts[3,3]
replace V_pr=bt in `x'
scalar bt = bts[4,4]
replace V_t_pr=bt in `x'
scalar bt = bts[5,5]
replace V_enpp=bt in `x'
scalar bt = bts[6,6]
replace V_t_enpp=bt in `x'
scalar bt = bts[7,7]
replace V_vote=bt in `x'
scalar bt = bts[8,8]
replace V_t_vote=bt in `x'
scalar bt = bts[9,9]
replace V_partyold=bt in `x'
scalar bt = bts[10,10]
replace V_t_partyold=bt in `x'
scalar bt = bts[11,11]
replace V_t_niche=bt in `x'
scalar bt = bts[12,12]
replace V_gov=bt in `x'
scalar bt = bts[13,13]
replace V_t_gov=bt in `x'
scalar bt = bts[14,14]
replace V_inc=bt in `x'
scalar bt = bts[15,15]
replace V_t_inc=bt in `x'
scalar bt = bts[16,16]
replace V_cons=bt in `x'
scalar bt = bts[17,17]
replace V_level3=bt in `x'
scalar bt = bts[18,18]
replace V_level2=bt in `x'
scalar bt = bts[19,19]
replace V_level1=bt in `x'

ereturn clear

}

drop if b_daysbeforeED==.
keep m M b_* se_* V_* _*_log_ae daysbeforeED enpp t_enpp pr t_pr vote t_vote partyold t_partyold t_niche gov t_gov inc t_inc cipollse cons

save MI_MLWin.dta, replace

** -> CALCULATE MI ESTIMATES BASED ON RUBIN'S RULES *****************************************************************************************************************

* See: https://bookdown.org/mwheymans/bookmi/rubins-rules.html

use MI_MLWin.dta, clear

gen level1=.
gen level2=.
gen level3=.

* SET NUMBER OF OBSERVATIONS (N)
gen n=359024

* SET NUMBER OF PARAMATERS (K)
gen k=19

foreach x of varlist daysbeforeED enpp t_enpp pr t_pr vote t_vote partyold t_partyold t_niche gov t_gov inc t_inc cipollse cons level1 level2 level3 {

* GENERATE POOLED PARAMETER ESTIMATE ('POOLED MEAN DIFFERENCE')
egen u_b_`x'=mean(b_`x')

* GENERATE WITHIN IMPUTATION VARIANCE
egen Vw_`x'=mean(V_`x')

* GENERATE BETWEEN IMPUTATION VARIANCE
* (i) DIFFERENCE BETWEEN EACH PARAMETER ESTIMATE AND THE POOLED MEAN DIFFERENCE
gen mdiff_`x'=(b_`x'-u_b_`x')^2

* (ii) SUM OF DIFFERENCES FROM (i)
egen sum_mdiff_`x'=total(mdiff_`x')

* (iii) AVERAGE OF DIFFERENCES FROM (i) & (ii)
gen Vb_`x'=sum_mdiff_`x'/(M-1)

* GENERATE TOTAL VARIANCE
gen Vt_`x'=Vw_`x'+Vb_`x'+(Vb_`x'/M)

* GENERATE POOLED STANDARD ERROR
gen u_se_`x'=sqrt(Vt_`x')

* GENERATE POOLED WALD TEST STATISTIC (ASSUMING PARAMETER VALUE UNDER THE NULL = 0)
gen Wald`x'=((u_b_`x'-0)^2)/(Vt_`x')

* GENERATE LAMBDA
gen lambda`x'=(Vb_`x'+(Vb_`x'/M))/Vt_`x'

* GENERATE RELATIVE INCREASE IN VARIANCE
gen riv`x'=(Vb_`x'+(Vb_`x'/M))/Vw_`x'

* GENERATE DF (OLD)
gen df_old`x'=(M-1)/((lambda`x')^2)

* GENERATE DF (OBSERVED)
gen df_obs`x'=((((n-k)+1))/(((n-k)+3)))*(n-k)*(1-lambda`x')

* GENERATE DF (ADJUSTED)
gen df_adj`x'=(df_old`x'*df_obs`x')/(df_old`x'+df_obs`x')

* DERIVE P-VALUE FOR POOLED MEAN DIFFERENCCE
gen absWald`x'=abs(Wald`x')
gen p_`x'=(2*ttail(df_adj`x',abs(Wald`x')))

}

keep m M b_* V_* u_* Wald* p_*

collapse (mean) u_* p_*

save MI_MLM_FullModel.dta, replace



